Implementación del Gradient Descent Algorithm para encontrar parámetros de un modelo de regresión logÃstica¶
Recordemos que un modelo de regresión logÃstica estamos buscando una combinación de variables que modele el log-odds de la siguiente forma,
$$ \text{log} \biggl{(}\frac{P_i}{1-P_i}\biggl{)} = \beta_0 + \sum_{j=1}^{p} \beta_j x_{i,j} $$Donde $p$ representa la cantidad de variables (columnas, caracterÃsticas) en los datos, y $i$ es cada una de las filas, tal que $i \in \{1,\dots,n\}$ con $n$ el tamaño de muestra.
En un modelo de regresión tradicional podemos utilizar métricas como el error cuadrático medio (MSE) o error absoluto medio (MAE) para encontrar los parámetros que minimizan dicha expresión. En el caso de la regresión logÃstica dichas métricas no sirven ya que nuestra variable respuesta tiene la forma $y_i \in \{0,1\}$ considerando un problema de clasificación binomial.
El modelo clásico de regresión logÃstica hace uso de la función logÃstica para modelar las clases / categorÃas de la variable respuesta, la cual tiene la siguiente forma,
$$ P_i(\beta, x_i) = \frac{1}{1+\text{exp}(-(\beta_0 + \sum_{j=1}^{p})\beta_j x_{i,j})} $$Para entender la función logÃstica, consideremos una forma sencilla de su expresión,
$$ f(z) = \frac{1}{1+\text{exp}(-z)} $$Si $z \rightarrow ∞$ entonces $f(z) = 1 $
Si $z → -∞$ entonces $f(z) = 0$
Si $z = 0$ entonces $f(z) = \frac{1}{2}$
Introduciendo una función de costo¶
Tal como se mencionó anteriormente, la función de costo asociada al problema de clasificación no puede ser del conjunto de expresiones utilizadas en regresión para respuesta continua, para este caso se introduce el concepto de entropia cruzada (cross-entropy), con esto podemos establecer una métrica que mide qué tanto se equivoca nuestro modelo en el problema de clasificación
$$ Costo = -\frac{1}{n} \sum_{i} y_i\text{log}(\hat{y}_i) + (1-y_i) \text{log}(1-\hat{y}_i) $$Donde $\hat{y}$ representa los valores predichos por el modelo, considerando:
$$ \hat{y}_i = \Phi(\beta_0 + \sum_{j=1}^{p} \beta_j x_{ij}) $$$$ \Phi(\hat{y}_i) = \frac{1}{1+\text{exp}(-\hat{y}_i)} $$Esto nos permite realizar una implementación del algoritmo de gradiente para la estimacion de parámetros. Para simplicidad del ejercicio vamos a trabajar con un dataframe con dos variables $x_1$ y $x_2$ y una respuesta $y \in \{0,1\}$.
AsÃ, las funciones de actualización quedan definidas por:
$$ \beta_j = \beta_j - \alpha \times \biggl{(}\frac{\partial Costo}{\partial \beta_j}\biggl{)}, \forall j \in \{1,2\} $$$$ \beta_0 = \beta_0 - \alpha \times \biggl{(} \frac{\partial Costo}{\partial \beta_0} \biggl{)} $$Sin entrar en el detalle matemático (el cual puede ser relativamente complejo ya que se ocupa derivada implicita, regla de la cadena, etc.), las derivadas parciales tienen la siguiente forma,
$$ \frac{\partial Costo}{\partial \beta_j} = (\hat{y} - y) \times x_j $$$$ \frac{\partial Costo}{\partial \beta_0} = (\hat{y} - y) $$Lo que nos permite obtener de forma más explicita las funciones de actualización
$$ \beta_1 = \beta_1 - \alpha\times (\hat{y} - y) \times x_1 $$$$ \beta_2 = \beta_2 - \alpha\times (\hat{y} - y) \times x_2 $$$$ \beta_0 = \beta_0 - \alpha\times (\hat{y} - y) $$Notar que esta implementación es generalizable para $p$ variables
Implementación en R¶
require(ggplot2)
require(dplyr)
# Simulamos datos para clasificación
# Número de filas en la simulación
n_datos <- 100
df <- data.frame(X1 = c(rnorm(n = n_datos/2, mean = 1, sd = 1),
rnorm(n = n_datos/2, mean = 5, sd = 1)),
X2 = c(rnorm(n = n_datos/2, mean = 1, sd = 1),
rnorm(n = n_datos/2, mean = 5, sd = 1))) %>%
# Agregamos una columna con 0 y 1 como variable respuesta
mutate(y = c(rep(0, n_datos/2), rep(1, n_datos/2)))
# Reordenamos los datos
df <- df[sample(1:nrow(df), size = nrow(df)), ]
# Opciones de visualizacion
options(repr.plot.width = 5, repr.plot.height = 5, repr.plot.res = 300)
# Visualizacion
df %>%
# Cambiamos a factor la respuesta
mutate(Label = as.factor(y)) %>%
ggplot(aes(x = X1, y = X2, col = Label)) +
geom_point(aes(shape = Label), size = 2.2) +
xlab('Variable 1') +
ylab('Variable 2') +
# Definimos colores manuales
scale_color_manual(values = c("deeppink","darkslategray4")) +
# Definicion de caracteristicas varias del grafico
theme(panel.background = element_rect(fill = 'beige', colour = 'black'), # Color de fondo
plot.background = element_rect(fill = 'beige'), # Color de fondo
legend.background = element_rect(fill = 'beige'), # Color de fondo leyenda
panel.grid.major = element_blank(), # Eliminamos la grilla
panel.grid.minor = element_blank(), # Eliminamos la grilla
axis.text.x = element_blank(), # Eliminamos texto de ejes X e Y
axis.text.y = element_blank(),
axis.ticks.x = element_blank(), # Eliminamos ticks de los ejes
axis.ticks.y = element_blank(),
legend.title = element_text(size = 15), # Tamaño del titulo del grafico
legend.text = element_text(size = 15), # Tamaño del texto del grafico
legend.key.size = unit(2,"line") # Tamaño de los simbolos en la leyenda
)
# Algoritmo
# Ratio de aprendizaje
alpha <- 0.001
# Epocas de entrenamiento
epocas <- 1000
# Definicion inicial de los parametros
beta = 0
beta_1 = 0
beta_2 = 0
# Vector para guardar los valores de la pendiente y el intercepto
slopes <- rep(NA, epocas)
ints <- rep(NA, epocas)
betas <- rep(NA, epocas)
betas_1 <- rep(NA, epocas)
betas_2 <- rep(NA, epocas)
for (i in 1:epocas){
logit <- beta + beta_1 * df$X1 + beta_2 * df$X2
y <- 1 / (1 + exp(-logit))
beta <- beta - sum(y - df$y) * alpha
beta_1 <- beta_1 - sum((y - df$y) * df$X1) * alpha
beta_2 <- beta_2 - sum((y - df$y) * df$X2) * alpha
slopes[i] <- -beta_1 / beta_2
ints[i] <- -beta / beta_2
betas[i] <- beta
betas_1[i] <- beta_1
betas_2[i] <- beta_2
}
En lo anterior estamos calculando los slopes e interceptos diviendo los coeficientes obtenidos, esto se debe a lo siguiente,
$$ 0 = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} $$Si queremos visualizar, tenemos que dejar una de las variables en función de la otra, para este caso despejamos $x_{2}$ en función de $x_{1}$
$$ x_{2} = -\frac{\beta_0}{\beta_2} - \frac{\beta_1}{\beta_2} x_{1} $$Lo que corresponde a una ecuación de la recta tradicional de la forma $y = mx + b$
betas_df <- data.frame(n = 1:epocas, beta = betas, beta_1 = betas_1, beta_2 = betas_2)
betas_df %>%
ggplot(aes(x = n)) +
geom_line(aes(y = beta), col = 'green') +
geom_line(aes(y = beta_1), col = 'red') +
geom_line(aes(y = beta_2), col = 'red') +
xlab('Epocas') +
ylab('Parámetros beta') +
annotate("text", x = tail(epocas)[1] * 0.95, y = tail(betas)[1] * 0.9, label = "beta 0") +
annotate("text", x = tail(epocas)[1] * 0.95, y = tail(betas_1)[1] * 1.1, label = "betas") +
theme(panel.background = element_rect(fill = 'beige', colour = 'black'), # Color de fondo
plot.background = element_rect(fill = 'beige'), # Color de fondo
legend.background = element_rect(fill = 'beige'), # Color de fondo leyenda
panel.grid.major = element_blank(), # Eliminamos la grilla
panel.grid.minor = element_blank(), # Eliminamos la grilla
legend.title = element_text(size = 15), # Tamaño del titulo del grafico
legend.text = element_text(size = 15), # Tamaño del texto del grafico
legend.key.size = unit(2,"line") # Tamaño de los simbolos en la leyenda
)
parametros <- data.frame(n = 1:epocas, slopes = slopes, ints = ints)
parametros %>%
ggplot(aes(x = n)) +
geom_line(aes(y = slopes), col = 'red') +
geom_line(aes(y = ints), col = 'blue') +
xlab('Epocas') +
ylab('Parametros') +
annotate("text", x = tail(epocas)[1] * 0.9, y = tail(slopes)[1] * 1.2, label = "Slope") +
annotate("text", x = tail(epocas)[1] * 0.9, y = tail(ints)[1] * 0.95, label = "Intercept") +
theme(panel.background = element_rect(fill = 'beige', colour = 'black'), # Color de fondo
plot.background = element_rect(fill = 'beige'), # Color de fondo
legend.background = element_rect(fill = 'beige'), # Color de fondo leyenda
panel.grid.major = element_blank(), # Eliminamos la grilla
panel.grid.minor = element_blank(), # Eliminamos la grilla
legend.title = element_text(size = 15), # Tamaño del titulo del grafico
legend.text = element_text(size = 15), # Tamaño del texto del grafico
legend.key.size = unit(2,"line") # Tamaño de los simbolos en la leyenda
)
# Ploteamos la regla de clasificacion
df %>%
mutate(Label = as.factor(y)) %>%
ggplot(aes(x = X1, y = X2, col = Label)) +
geom_point(aes(shape = Label), size = 2.2) +
geom_abline(aes(slope = tail(slopes)[1], intercept = tail(ints)[1]), color = 'red', linetype = 'dashed', size = 1.1) +
xlab('Variable 1') +
ylab('Variable 2') +
scale_color_manual(values = c("deeppink","darkslategray4")) +
theme(panel.background = element_rect(fill = 'beige', colour = 'black'),
plot.background = element_rect(fill = 'beige'),
legend.background = element_rect(fill = 'beige'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank())
Generación del GIF¶
Para generar el GIF vamos a tomar un gráfico para distintos valores de los interceptos y pendientes calculados, luego los juntaremos mediante gifski y visualizaremos
# setwd('Acá deberÃa ir tu working directory')
epocas_gif <- seq(from =1, to = epocas, by = 5)
for (i in epocas_gif){
gg <- df %>%
mutate(Label = as.factor(y)) %>%
ggplot(aes(x = X1, y = X2, col = Label)) +
geom_point(aes(shape = Label), size = 2.2) +
geom_abline(aes(slope = slopes[i], intercept = ints[i]), color = 'red', linetype = 'dashed', size = 1.1) +
xlab('Variable 1') +
ylab('Variable 2') +
scale_color_manual(values = c("deeppink","darkslategray4")) +
theme(panel.background = element_rect(fill = 'beige', colour = 'black'),
plot.background = element_rect(fill = 'beige'),
legend.background = element_rect(fill = 'beige'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank())
ggsave(sprintf("myplot%03d.png", i), gg)
}
# install.packages(gifski)
require(gifski)
png_files <- Sys.glob("myplot*.png")
gifski(
png_files,
'myanimation.gif',
width = 400, height = 400,
delay = 1/10
)
# Eliminamos los archivos png guardados para dejar solo el GIF animado
file.remove(png_files)
Lo anterior entrega un GIF de la siguiente forma, notar cómo se mueve el intercepto y la pendiente de la recta que realiza el corte entre las categorÃas
Lo anterior puede ser extendido a más variables, para mantener la visualización añadamos una nueva variable
# ----------------------------------
# Tres dimensiones -----------------
# ----------------------------------
# Agregamos una nueva variable
df <- df %>% mutate(X3 = c(rnorm(n = n_datos/2, mean = 1, sd = 1),
rnorm(n = n_datos/2, mean = 5, sd = 1)))
require(scatterplot3d)
alpha <- 0.001
# Epocas de entrenamiento
epocas <- 1000
# Definicion inicial de los parametros
beta = 0
beta_1 = 0
beta_2 = 0
beta_3 = 0 # Nuevo parametro dado que tenemos una nueva variable
# Vector para guardar los valores de la pendiente y el intercepto
slopes_1 <- rep(NA, epocas)
slopes_2 <- rep(NA, epocas) # un slope extra
ints <- rep(NA, epocas)
for (i in 1:epocas){
logit <- beta + beta_1 * df$X1 + beta_2 * df$X2 + beta_3 * df$X3
y <- 1 / (1 + exp(-logit))
beta <- beta - sum(y - df$y) * alpha
beta_1 <- beta_1 - sum((y - df$y) * df$X1) * alpha
beta_2 <- beta_2 - sum((y - df$y) * df$X2) * alpha
beta_3 <- beta_3 - sum((y - df$y) * df$X3) * alpha # Calculo del tercer parametro con actualizacion
slopes_1[i] <- -beta_1 / beta_3
slopes_2[i] <- -beta_2 / beta_3
ints[i] <- -beta / beta_3
}
aux <- 0
for (i in epocas_gif){
# Dado que vamos a plotear con librerias distintas a ggplot tenemos que hacer ajustes en los graficos
graphics.off() # Dejamos el graphics.off() para guardar el plot
png(filename = sprintf("myplot%03d.png", i), width = 8, height = 6, units = "in", res = 300) # Nombre de los archivos a generar
colors <- c("deeppink","darkslategray4") # Colores para el 3D plot
colors <- colors[as.factor(df$y)] # Colores como vector
s3d <- scatterplot3d(df[,c('X1','X2','X3')], # Scatterplot3D con las columnas a plotear
color = colors,
pch = 16,
box = FALSE,
# Limites customizados
xlim = c(-10,10),
ylim = c(-10,10),
zlim = c(-10,10),
# Angulo que gira dinamicamente
angle = 310+aux/5)
# Agregamos el plano con el intercepto y coeficientes
s3d$plane3d(Intercept = ints[i],
x.coef = slopes_1[i],
y.coef = slopes_2[i], lty = 'dashed', lty.box = 'solid')
dev.off()
aux <- aux + 1
}
png_files <- Sys.glob("myplot*.png")
gifski(
png_files,
'myanimation_3D.gif',
width = 400, height = 400,
delay = 1/10
)
file.remove(png_files)
Lo anterior genera el GIF mostrando cómo se ajustan los parámetros de la regresión logÃstica considerando las tres variables $x_1, x_2, x_3$. De forma similar al plot en 3D tenemos que ajustar el plano generado mediante lo siguiente,
$$ x_3 = -\frac{\beta_0}{\beta_3} -\frac{\beta_1}{\beta_3} x_1 - \frac{\beta_2}{\beta_3}x_3 $$